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Abstract 

We present preliminary results on a novel numerical method describing wave 
propagation in slowly non-uniform media. Following Huygens-Fresnels prin¬ 
ciple, we model the wavefront as an array of point sources that emit wavelets, 
which interfere. We then identify a set of new points where the electric field 
has equal phase. In fact, without losing generality, we find zeros of the elec¬ 
tric held, by means of the bisection method. This obviously corresponds to 
a specihc phase-advance, but is easily generalized, e.g. by phase-shifting all 
sources. The points found form the new wavefront. One of the advantages 
of the method is that it includes diffraction. Two examples provided are 
diffraction around an obstacle and the hnite waist of a focused Gaussian 
beam. Refraction is also successfully modeled, both in slowly-varying media 
as well as in the presence of discontinuities. The calculations were performed 
in two dimensions, but can be easily extended to three dimensions. We also 
discuss the extension to anisotropic, birefringent, absorbing media. 

Keywords: Huygens-Fresnel principle, wavefront tracing, ray tracing, full 
wave, computational electromagnetism, diffraction 
PACS: 42.15.Dp, 42.25.Bs, 42.25.Fx, 42.25.Gy, 42.25.Lc 


1. Introduction 

Huygens stated in 1690 that “each element of a wavefront may be re¬ 
garded as the centre of a secondary disturbance which gives rise to spherical 
wavelets”. He added that “the position of the wavefront at any later time is 
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the envelope of all such wavelets” BE]- This result is referred to as Huy¬ 
gens’ theorem or construction. Indeed we can use it to construct subsequent 
wavefronts, if we ascertain some subtleties in what “envelope” means in this 
context. Such subtleties were eventually addressed by Fresnel lasiE], as dis¬ 
cussed below. Also, the construction should be regarded as a mathematical 
abstraction that correctly reproduces the physics without necessarily being 
physically rigorous. Treating the wavefront elements as actual centers of sec¬ 
ondary disturbances might be appropriate in some physics areas, but should 
not be taken literally in optics: obviously “light does not emit light; only 
accelerating charges emit light” [B], and the construction should be put in 
the context of what was known and understood in 1690. Yet, a principle of 
equivalence guarantees the validity of Huygens’ results: given two adjacent 
domains, 1 and 2, and given a wave in domain 1 approaching domain 2, it 
is possible to dehne conditions at the boundary between 1 and 2 (charges 
oscillating at the proper frequency, amplitude and direction) such that the 
wave excited in domain 2 by these boundary sources is indistinguishable from 
the wave entering from domain 1. Note that the charges in question oscillate 
as if responding to the wave in domain 1. The only imperfection in Huygens 
construction is that it pictured point-emitters, whereas elementary dipoles 
are more appropriate, as discussed in Ref. [7] and in Secj^ 

In 1816 Fresnel postulated that Huygens’ secondary wavelets interfere 
with each other 00]. At this point, we can look at the interference pattern, 
and identify the next wavefront as the locus of points with the same phase. 

The combination of Huygens’ construction with the principle of inter¬ 
ference is called Huygens-Fresnel principle [2] or, often, simply Huygens’ 
principle, and is an important principle in the theory of diffraction. 

Hence, not surprisingly, the principle has been used for didactic physics 
widgets illustrating diffraction around an obstacle, diffraction through one 
or more slits [8], and the diffraction of a Gaussian beam (namely, its hnite 
waist). The principle also accounts for refraction, provided wavelets expand 
with different velocities in different media, in inverse proportion with the 
respective indices of refraction. For this reason, Huygens’ principle has also 
been applied to Snell’s law |9]. 

Recently Huygens’ principle was also invoked in the design of metamate¬ 
rials with new beam shaping, steering and focusing capabilities |10j . 

In addition, the Huygens’ construction and Huygens-Fresnel’s principle 
inspired numerical methods, briefly reviewed in the next two subsections. 
This is also not surprising, on considering that Huygens’ original drawings 
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PP basically outlined an iterative graphical solver, vaguely reminescent of a 
modern iterative solver on a computer. 

1.1. Seismic wavefront tracing methods 

Huygens’ principle was considered in Ref. m, but discarded due to the 
lack of physics rigor mentioned above [ 6 ]. 

Wavefronts are “traced” or “tracked” in seismology on the basis of ray¬ 
tracing results or of geometrical-optics-approximated calculations of shortest 
travel-time to points on a grid [12], which in optics is equivalent to Fermat’s 
principle of least time. In particular, among grid-based methods, fast march¬ 
ing methods [13] and some earlier works [S] are close to Huygens’ idea, in 
that they treat the wavefront as a moving interface. Ray tracing ignores 
diffraction by dehnition, but calculations of shortest travel-time interrogate 
all points in a domain, including points around an obstacle. This might look 
promising for the inclusion of diffraction effects. Yet, it should be pointed out 
that these methods typically solve the Wentzel-Kramers-Brillouin (WKB) ap¬ 
proximated eikonal equation. Therefore, the wave field is only accurate in 
the high-frequency, short-wavelength limit. 

An interesting and elegant method is developed in Refs. [T5l [T 6 ] . Again, 
the starting point is the eikonal equation, but manipulated into the expres¬ 
sion: 

[x - x( 7 , + [y- 2/(7, 0)]^ + [^ - ^(7, 0)]^ = 0), (1) 

where 7 and 0 are curvilinear coordinates on the wavefront, x{'y,(j)), 1 /( 7 , 0 ), 
2 :( 7 , 0) are the Cartesian coordinates of a (source) point on the wavefront, 
and x,y,z the coordinates of another point, which can be pictured as a re¬ 
ceiver. Eqj^ describes a sphere centered at a:( 7 , 0 ), 1 /( 7 , 0 ), ^( 7 , 0 ), of radius 
r( 7 , 0 ) = u( 7 , 0 )r, where u is a velocity. As time r progresses, the sphere 
expands. Differentiating Eqjl] eventually leads to an explicit hnite-difference 
scheme: the coordinates of a point on the new wavefront are functions of 
the local r and of the coordinates of some points on the current wavefront. 
In particular, stencils of three and hve points are used respectively in the 
2D and 3D problem. The analogies with Huygens’ construction are evident. 
However, Eqjl] contains no information on the strength of individual point- 
sources. Also, each new point is only informed by three or five points in 
its vicinity, whereas diffraction integrals imply that every point on the new 
wavefront depends on every point on the old one. Hence, although the nu¬ 
merical technique is highly stable [T51ITB] . small phase-increments might be 
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required in between consecutive wavefronts, in order for the restriction to 
three or hve points to be a good approximation. 


1.2. Other wavefront tracing and Huygens-related methods 

Transmission Line Matrix Modeling ini maps the problem of wave prop¬ 
agation in a medium into the problem of electrical signal propagation in a 
3D array of shunt and series nodes. It can be shown that this is equivalent 
to Huygens’ mechanical model of light as a perturbation of the Ether, origi¬ 
nated by intense heat and transmitted by elastic shocks from one particle to 
the next, in all directions im. It is also equivalent to a many-body scatter¬ 
ing problem, in which each node scatters light emitted by other nodes. It is 
possible to include diffraction by proper formulation of the scattered wave, 
or by adding ad hoc dielectrics to the array mentioned above [18] . 

Face offsetting [19] is a Lagrangian algorithm that treats an interface as 
a surface-mesh, propagates the mesh faces and, from their envelope, recon¬ 
structs the next surface, and its vertices. The method is based on a gener¬ 
alization of Huygens’ original construction and envelope idea, but without 
Fresnel’s addition of wavelet-interference. 

Diffraction by rough surfaces (thus, ultimately, Huygens-Fresnel princi¬ 
ple) is the physical explanation for iridescent colors. However, their render¬ 
ing in computer graphics is simply obtained by angle-dependent coloration 
of surfaces in ray tracings m- 

Recent algorithms 1211 utilize the Huygens-Kirchhoff integral to integrate 
many locally valid asymptotic Green functions into a globally valid asymp¬ 
totic Green function for the Helmoltz equation in inhomogeneous media. 
However, they are formulated under geometrical optics approximations. 

Additional methods and potential numerical applications of Huygens’ 
principle are discussed in Ref. [22]. 

Note that “wavefront tracing” here refers to reconstructing several wave- 
fronts in a propagation problem, whereas “wavefront reconstruction” in large 
telescopes denotes the inverse problem (phase retrieval) of reconstructing, 
from sensor measurements, a single incoming wavefront, in order to decide 
how to make it planar by means of adaptive optics [23|. Yet, there might be 
synergies with the method proposed here, one step of which is indeed phase 
retrieval (Sec,2.3). Wavefront reconstruction includes diffraction. 
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1.3. Present work 

In summary, although several numerical works use Huygens construction 
of emitters and wavelets, few of them utilize the Huygens-Fresnel principle 
in its entirety, including wavelet interference (that is, diffraction). 

A very simple numerical method is presented in Secj^ of the present ar¬ 
ticle, that retains more than 3 or 5 points in the diffraction integrals, takes 
into account their different intensities and utilizes a zero-search to localize 
the points forming the new wavefront. Secj^ presents numerical examples 
obtained in a 2D isotropic non-uniform medium. 

The original motivation for this work stemmed from the search for wave 
propagation methods in magnetized plasmas, intermediate in physics content 
and computational cost between ray-tracing pTl |25l |26] and beam-tracing or 
paraxial Wentzel-Kramers-Brillouin (WKB) methods [?n |2H] on one hand, 
and full-wave solvers |2Hl ED] on the other. In turn, such need was prompted 
by unexpected results from the first full-wave study in the electron cyclotron 
frequency range [3T] -a high-frequency range, typically well-modeled by ray 
tracings. However, the method proposed here is not restricted to magnetized 
plasmas, and could be extended to generic 3D anisotropic birefringent media, 
as discussed in SecJH 


2. Numerical Method 
2.1. Governing equation 

Consider a wave in a medium, of wavenumber k = 271/X, where A is the 
wavelength in the medium. In general, k and A differ from the wavenumber 
and wavelength in vacuum. 

The Huygens-Fresnel formula for the complex field amplitude E in a point 
x' can be written as follows |32j : 


E(x') 



ik{l -I- cosO) 


—E(x)d^, 
r 


( 2 ) 


where x is the generic point on surface S', representing the initial wavefront. 
9 denotes the angle between the local wavefront-normal and the vector x'-x, 
of norm r. 

Strictly speaking, Eqj^is only correct for uniform media, which, however, 
can often also be treated analytically. Nonetheless, it is not unusual and 
perhaps more useful and interesting to also use Eqj^ [3S1 EH EH EH EZj or 
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other Huygens-based expressions [121 [T51 UHl Ell ES] in weakly non-uniform 
media (in the sense that these media can be treated as locally uniform on 
the lengthscale of a single numerical step, but could have different refractive 
properties in the following step, and the wave have a different k and A as a 
result). 

Said otherwise, the formula is still applicable if the steplength is prop¬ 
erly chosen, so that x' is “not too far” from any x, in the sense that L = 
(A^) ::§> |x'-x|. Here L is the lengthscale over which the medium (hence, 

the wavenumber k) varies. Note that the constraint on x' not being too far 
depends on L. Also note that the approximation does not necessarily depend 
on how L compares with the wavelength A=27r/fc. In other words, even if 
the medium varies signihcantly over a wavelength A, Eqj^can still be used 
to calculate E in a point x' close enough to S (“close enough” in the sense 
that the medium is practically uniform over that short distance). In fact, 
any ordering of the wavelength A with respect to L and |x'-x| is acceptable: 
greater than both, smaller than both, or intermediate. 

In 2D, the surface integral is replaced by a line integral, and the denom¬ 
inator by a/t, because the intensity of a cylindrical wave decays like 1 /r 



E(x') = 


1 -|- cos 6 


2A 


E(x)d/. 


(3) 


The integrand is more complicated than expected from Huygens’ simplistic 
point sources. To begin with, the elementary sources emit with different 
intensities in different directions. This was conjectured by Fresnel to recon¬ 
cile theory and experiment, and was rigorously calculated by Stokes [2]. It 
responds to the intuition that, given a source and two observers, one with 
a direct line of view, the other located around a corner, the former should 
receive more light. 

Another correction is the cosO/r term in Eqj^ This is needed to make 
Huygens-Formula valid also in the near-held, and consistent with the more 
general and rigorous Kirchhoff integral theorem [7]. Its physical meaning is 
that the elementary sources should actually be radiating dipoles. 

Indeed, dipoles are routinely used in modeling meta-lenses, frequency- 
selective surfaces and other diffraction-based devices gni 011112]. 

Note that Eqsj^ and are not equations to solve. They are rather solu¬ 
tions, in 3D and 2D, of the Helmholtz equation governing electromagnetic 
wave propagation in uniform (and, with some approximation, slowly non- 
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uniform) media. With special care, they can also describe propagation in 
strongly non-uniform media (see Figs|^and|^and accompanying discussion). 

These solutions are constructed as superpositions of point-to-point solu¬ 
tions, propagating from a point-emitter to a point-observer. The treatment 
presented in the present paper permits to efficiently calculate integrals or 
ID in realistic domains with obstacles and non-uniformities. Its outcome are 
wavefronts reconstructed with high precision -strongly sub-wavelength. Yet, 
these hnely reconstructed wavefronts can be sought at large distances from 
each other (large compared with the wavelength, but small compared with 
the lengthscale of inhomogeneity). 

The treatment presented is easily generalized to other problems, provided 
the point-to-point propagator (Green’s function) is known, and the superpo¬ 
sition principle is valid. 

2.2. Wavefront construction by zero localization 

Given a wavefront S', our objective is to compute the integral in EqJ^to 
identify a set of points x.'ab (with a=l,... m and b=l,.. .n) where E has equal 
phase. Without losing generality, we can search for 3fJ(E)=0: if E = Ae*^, 
then 3fJ(E)=0 implies either zero amplitude (A = 0) or a specihc phase 
(cos 0=0). However, the trivial solution of zero amplitude is easily rejected: 
zero amplitude means no wave, and if there is no wave in x', presumably 
there is no wave in a neighborhood of x'. By contrast, 3fJ(E)=0 as due to 
cos 0=0 can only be true on a particular wavefront. 

This search for a specihc phase (0 = 7r/2 modulo vr) is easily generalized 
to a different phase by shifting all sources in the integrand by a certain 50. 
The points found form the new wavefront. 

Equivalently, one could search for the locus of points where the phase of 
the total electric held, arctan[9^(E)/jR(E)] equals the phase of interest, 0. 
This can also be casted as a zero-search problem. 

The existence of multiple solutions is prevented by restricting the search 
to intervals of length A/2, where A is the wavelength in the medium. It should 
be clarihed that the length of the search-interval has nothing to do with its 
positioning-, the search-interval can be placed close to the wavefront S or far 
from it, depending on the desired resolution (in the sense of inter-wavefront 
spacing) and on the validity of Eqj^ at that distance from S (in other words, 
is it still |x' — x| -C [|^(x)] ^ for any x' on the A/2 long search-interval and 
for every x on surface S7). Note that the wavefront-to-wavefront distance 
depends on the positioning of the search-intervals, not on their length. An 
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important consequence is that acceptable wavefront-to-wavefront distances 
(thus, computational costs) do not depend on A; they only depend on the 
inhomogeneity of the medium. 

Depending on where the user places the search-intervals, the algorithm 
will trace the wavefront at phase 0 = 7r/2, or 37r/2, or 57r/2 etc. from the 
current one. As discussed, this is because sources have zero initial phase in 
Eq§ Nonetheless, different wavefront-to-wavefront phase-increments (dif¬ 
ferent from odd multiples of vr/2) can also be obtained. For this purpose 
we can perturb the phase at the sources (e*^^E(x) —>■ e*^'’’''®'^E(x)) while still 
searching for zeros at the receivers (3fJ[E(x')] = 0). 

At least two approaches can now be envisioned for either zero search (that 
is, of 3?(E) or arctan[A(E)/3f?(E)] — 0). 

A possible approach is to locally solve Maxwell’s equations in such vol¬ 
ume, in presence of point-sources located on the current wavefront. Com¬ 
putational electromagnetics methods, including but not limited to hnite dif¬ 
ferences and hnite elements, could repetitively solve the “local” problem in 
front of the present wavefront. However, in general this approach offers no 
major advantage compared to solving the problem once in the “global” do¬ 
main. The approach is only advantageous if the union of the local domains 
is smaller than the global domain. 

More interestingly, and more in line with the spirit of the Huygens-Fresnel 
principle, we can let the original point-sources interfere and we can localize 
iso-phase points to construct the new wavefront. In particular, we can look 
for loci of destructive interference (zeros). 

2.3. Varying the phase-step 

The zero-search is conducted on a A/2 long segment. Depending on the 
positioning of the search intervals, different zeros could be found, correspond¬ 
ing to a different wavefront or, which is the same, to a different phase-advance 
A0 (modulo tt) with respect to the previous front. It is necessary to con¬ 
sistently place the search-intervals in such a way that they bracket zeros all 
belonging to the same wavefront. In other words, the search intervals need 
to “bracket” the unknown wavefront with a precision of A/2. 

This provides a simple knob to vary the phase-advance from one wavefront 
to the next: it is sufficient to consistently move all search-intervals farther 
from or closer to the current wavefront to vary the phase-advance by utt, 
where n is an integer. 


For finer adjustments of 0, one can use complex E in the integrand of 
Eqj^ or, equivalently, multiply the integrand by which is equivalent to 
phase-shifting the sources. 

2.4- Amplitude calculation and initialization of following iteration 

If the wavefront was constructed as a locus of zeros (3fJ(E)=0), it would 
be unclear what the wave amplitude E is in each point. This information 
is needed because the hnal wavefront of a given step is the initial wavefront 
for the following one, and its points need to be properly initialized with the 
correct amplitude. 

For this, it will be necessary to either compute Q'(E) in the wavefront 
points (if 3fJ(E) is already known), or re-compute 3fJ(E) in those same points, 
but with all sources phase-shifted by 0 7 ^ utt. Note that the wavefront points 
have already been identihed and do not need to be identihed again. 

Also note that this helps discarding “false zeros” of the wave-held. By 
this we mean that, in addition to zeros (3^(E) ~0) of the oscillating, non- 
negligible wave-held (|E| 7 ^ 0 ), there could be trivial zeros where 3fJ(E) ~0 
and |E| ~0. Such zeros are located away from the wave beam, and are not 
relevant to wavefront identihcation. 

2.5. Algorithm 

The algorithm -in 2D to hx the ideas, in Cartesian coordinates x and y, 
but easily generalized to 3D- can be recapitulated and illustrated as follows. 

Step 1 : Discretize the current wavefront as a ID array of points 
(a=l,...m) emitting diherent intensities E^ (color-coded in Figj^), dis¬ 
tributed according to the beam-pattern. 

The very hrst array can be at the antenna, phased array or last mirror 
and can have complicated shape or phase-pattern. 

Step 2 ; Dehne A/2 long search-intervals for the points x'^ (A=l,. .. m) 
that will form the next wavefront. 

For simplicity they can have the same coordinates y as the previous points 
Xa, and only span x (Figj^). In fact, they can be obtained by displacing all 
points Xa by the same amounts Ax and Ax -|- A/ 2 . Avoiding to place such 
intervals too closely (Ax A) to the original wavefront will allow discarding 
the dipole correction. 

Step 3 : For A going from 1 to m: 

a) Consider the search-interval for x'^. 
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b) Evaluate the field in a point on that interval, where c denotes 
an iteration index in the zero-search, e.g. by bisection |13]. Note 
that E^c (of fixed A) is obtained by interference of all point-sources 
(a=l,... m) in FigjJ:, each one of amplitude E^, at distance „ from 
the target and inclination 9Aa,c- 



This is a discretization of Eqj^ In its present form, the number of field 
calculations in m observers due to m emitters scales like However, 
unequally spaced Fast Fourier Transform methods |11] could make it 
scale like mlogm. Additionally, Eqj^ can be restricted to values of a 
(typically in the neighborhood of A) yielding non-negligible contribu¬ 
tions (dashed lines in Figj^). 

c) Repeat step b) for various x'^^ “suggested” by the zero-finder in its con¬ 
vergence to the zero. Repeat until |3fJ(E^c)| < e, where e is a prescribed 
tolerance. 

d) Set point x'^ equal to last x'^^. 

e) Set field-amplitude in x'^ equal to last A(Eac) (Figj^). 

Step 4: Repeat steps 1-3 until absorption, or until exiting from com¬ 
putational domain, or until other criterion is met. Use the final points and 
amplitudes of step 3 as initial points and amplitudes for step 1. 

The 2D algorithm is easily generalized to 3D by adding a subscript ^ and 
summing and looping over it. Eqj^ needs to be replaced by: 



3. Numerical examples 

The first test consisted of aiming a plane wave in the x direction, with 
periodic boundary conditions in ?/, effectively simulating an infinitely wide 
planar wavefront. This very simple test worked as expected: the wavefront 
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Figure 1: Illustration of (a) step 1, (b) step 2 and (c-d) step 3 of the algorithm described 
in Sec 2.5 for tracing a new wavefront based on the current one. 


propagated in the x direction and remained planar (Figj^). This is in obvi¬ 
ous agreement with the trivial analytical solution. 

Next, an obstacle was placed at the bottom left of the computational box 
(Figj^). No periodic boundary conditions were used here, but additional 
emitters (“ghost cells”) were added at large y, outside of the computational 
box of interest. These acted as point-sources, which were taken into account 
in localizing new zeros and constructing new wavefronts, as well as calculating 
the amplitudes in wavefront points, for proper initialization of each step. 
Figj^clearly exhibits diffraction around the obstacle, as expected. At shallow 
angles past the obstacle, as expected, the wave held is weak and its zeros are 
not shown. 

The next problem modeled was refraction in a non-uniform medium whose 
index of refraction varies slowly with space (ten wavelengths, in the case of 
Figj^. The wavefront is planar, oblique. It extends outside of the computa¬ 
tional domain by means of ghost cells on its top and bottom. As expected, 
the wavefront changes its direction of propagation according to the index of 
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Figure 2: (a) Simple test of tracing planar wavefronts in uniform medium, in the absence 
of obstacles, (b) Initially planar wavefronts diffract around obstacle as expected. 
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Figure 3: Wavefronts refracting in a medium of (a) refractive index N varying over ten 
wavelengths or (b) discontinuous. 


refraction, in agreement with the continous limit of Snell’s law. The wave¬ 
length also contracts or expands accordingly, as it can be recognized from 
the wavefront spacing. 

Discontinuities in the refractive index require special care. The reason 
is that the generic wavefront crossing a discontinuity has some points in 
medium 1, and some in medium 2. If E is being evaluated in a point in 
medium 2, the contributions from points in medium 1 will travel initially 
at speed c/Ni and then, once reached medium 2, at speed C/N 2 , where Ni 
and N 2 are the refractive indices in the two media. Due to the principle 
of least time, the information does not travel on a straight line, but on two 
consecutive segments of different inclinations, obeying Snell’s law. Therefore, 
instead of a single Green’s function relating emitter i with observation point 
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Figure 4: Illustration of why, in presence of discontinuities, a single Green’s function 
cannot propagate the information from an emitter in medium 1 to a receiver in medium 
2, and an intermediate point is needed at the boundary, determined by Snell’s law. 


j, 

ik-(xj-Xi) 

7 - 1 ’ ^ 6 ) 

|Xj -Xi| 

we should introduce an intermediate point x^ at the boundary, and use two 
propagators, from Xi j in medium 1 to x^ jj, and then from x;, to X 2 j in 
medium 2: 

„ik2-(x2,j-X6,ij) piki-(xi,,ij-xi,i) 

1—=—ri—=—r 

1^2,j- 1^6,^1,^1 

For a given pair of xij and X 2 j, ^bij is uniquely determined by Snell’s law 
(Figg. 

The wavefront bends as expected from Snell’s law (Figj^), in agreement 
with the analytical result. Corrugations similar to Figj^, but stronger, are 
ascribed to a lack of ghost cells, and will be the subject of future work. The 
construction of search-intervals will also be improved, by better bracketing 
of the zeros of interest. 

Related to refraction by a single discontinuity is refraction by two con¬ 
secutive ones. Wavefronts experience this when crossing a lens, and are 
observed to deform as expected, for example in the case of a convergent lens 
of refractive index iV=1.2 (Figj^. 

Propagation through a single lens can be easily generalized to propagation 
through a random medium consisting of variously shaped “blobs” or grains 
immersed in a background material of different refractive index. 

Finally, a Gaussian beam was simulated. Point-emitters were initialized 
on an initial curved wavefront, with intensities distributed according to a 
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Figure 5: Wavefronts focused by a convergent lens of refractive index N=\.2. 


Gaussian profile. Excellent results were obtained for a moderately focused 
Gaussian beam (Figj^), in good agreement with Gaussian optics theory 
(curves in Figj^) [15]. The geometrical optics solution is also shown and, as 
expected, is a good approximation away from the waist. 

Figj^does not show wave intensities smaller than 10“®, when normalized 
to intensities on the axis of the Gaussian beam. This corresponds to fractional 
electric helds of 10“^. Similar criteria were adopted in truncating the very 
edge of the wavefront diffracting around the obstacle in Figj^ 

At even smaller helds, the zero-search can fail and misplace portions of 
the wavefront. Although this error originates in typically uninteresting low-F 
regions, it should be noted that each wavefront is determined by the previous 
one. Hence, imprecisely traced portions of a wavefront might cause portions 
of “later” wavefronts to also be misplaced, or distorted. The ultimate result 
is that these errors can eventually affect hnite-F regions. This is typically not 
the case for moderately focused beams (Figj^), but can be the case for the 
strongly focused ones (Figj^), which will be the subject of future work. The 
issue could be solved by an improved, adaptive dehnition of the interval for 
zero-search by the bisection method. Additionally, here the zero search was 
performed for simplicity in the x direction (see search-intervals in Figj^). 
However, more customised search-intervals (orthogonal to the wavefront, or, 
better, parallel to group velocity) would aim at regions of high E and avoid 
low E and associated issues. 
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Figure 6: Wavefronts for (a) moderately and (b) strongly focused Gaussian beam, colored 
according to intensity normalized to on-axis value. One in every ten wavefronts is bold, 
as a guide for eyes. Also shown are the exact solutions for the contours at 1/e^ of central 
intensity (black curves on red background). 


4. Discussion 

By virtue of Fresnel’s principle, the method naturally includes diffraction. 
This is an advantage over ray tracings and truncated paraxial expansions used 
in beam-tracings. 

Another advantage of the method is that it only requires sub-wavelength 
resolution along the wavefront. Across the wavefront, however, the search 
for the next wavefront only requires a handful of points: if we assume the 
next wavefront to be one wavelength A away from the current one, then only 
log 2 100 iterations, -that is, 6 or 7- are needed to localize a point on the 
next wavefront with A/100 accuracy. This is an important advantage over 
methods needing 100 points or elements per wavelength. Moreover, the next 
wavefront can be searched at more than a wavelength from the previous one, 
resulting in a further speed-up. This is possible because complete diffraction 
formulas are used here, instead of near-held approximations. 

Going in more details, let us estimate the number of point-to-point calcu¬ 
lations Nops needed to trace A’wt wavefronts discretized by Ne points. Here 
by point-to-point we mean the contribution of a specihc point emitter to a 
specihc observation point (the argument of the sums in Eqj^. Localizing 
each point requires calculating the total held according to Eq|^ in Nbis trial- 
points. These are the points needed by the bisection method to converge to 
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a zero. In turn, each total field is calculated as the sum over a total of Ne 
emitters, in the worst case (i.e., on the net of truncations). Hence, a total of 
Nops = NeNwtNus point-to-point calculations is needed. 

Wavefront-spacing is not constrained by the wavelength, but needs to be 
much smaller than the non-uniformity length-scale normally to the wavefront, 
L_i_. Nevertheless, it is reasonable to assume that wavefront-spacing is a 
fraction / of the wavelength. Therefore, Nwt = -^a±//, where Nx± is the 
number of wavelengths in the direction orthogonal to the wavefront. In 
that case, reaching the orthogonal precision of a full-wave solver deploying 
Nppw points per unit wavelength imposes Wis = log 2 (/-^pp«))- In conclusion, 
Nops = NlNx^\og^{fN^p^)/f. 

As mentioned in Sec, 2.5[ unequally spaced Fast Fourier Transform meth¬ 
ods could further speed up the method by replacing Ne —)■ NE^ogNE- 


5. Future work 

In the original Huygens-Fresnel integral, k was a fixed scalar, but non¬ 
uniformity is easily taken into account by letting k = k{x), as done in Secj^ 

Also, the term in Huygens integrand assumed isotropy: A: is a scalar, 
meaning that wavelets expand in all directions x' — x with the same wave¬ 
length and same phase velocity u/k. 

In anisotropic media, however, global wavefronts as well as local wavelets 
expand at different velocities in different directions. Spherical wavelets are 
thus replaced by ellipsoids, elongated or compressed in the direction of anisotropy. 
Such ellipsoids do not preserve their aspect ratio, because some axes expand 
faster than others. It is easy to account for this by letting k = k{'x,x.' — x). 
The wavenumber now depends on the line-of-sight x' — x connecting the 
emitter with the observer (for example it depends on whether such direction 
is parallel, perpendicular or oblique to the direction of anisotropy), k is pro¬ 
vided by the dispersion relation for the wave of interest, in the polarization 
(mode) of interest, and in the direction x' — x of interest. An earlier treat¬ 
ment of anisotropy in terms of TE and TM modes |16] arrived to the same 
ellipsoidal wavelet representation. 

Finally, further generalizing k = k(x, x' — x, mode) to be a complex, 
mode-dependent vector, not necessarily parallel to x' — x, accounts respec¬ 
tively for absorption, birefringence and decoupling between phase and group 
velocity. The phase velocity is parallel to k, the group velocity describes 
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energy propagation from an array of points x to an array of points x'. Inci¬ 
dentally, phase and group velocity are decoupled also in ray tracings, allowing 
for non-parallel increments of ray position dr and of wavevector, dk. 

The underlying theory for k = k(x, x' — x, mode) (dispersion relation) is 
very well-established, but was never inserted in Huygens formulas. 

Further improvements can be realized by modifying the integrand in 
Eqsj2]j^by way of special functions: 

• Ref. [12] proposed to treat caustics by strategically ignoring some points 
of the old wavefront. In Ref. [33] it was proposed to address caustics by 
replacing with Airy functions the field radiated by point-sources. 

• Bessel functions, on the other hand, can model the diffraction of evanes¬ 
cent waves in the near field [32] . 

• Airy and parabolic cylinder functions can replace oscillating helds in 
some evanescent mode-conversion regions mi- 

Summary and Conclusions 

We explored the applicability of the Huygens-Fresnel principle to a new 
method for wave propagation in non-uniform media. 

Very briefly, the wavefront is discretized in an array of point-emitters 
of wavelets. A zero-search identifies observation points where the wavelets 
interfere destructively. Those zeros of held form an iso-phase surface, that is 
easily generalized to other phases. Then the process is repeated. 

Some very simple cases were successfully modeled in 2D, including diffrac¬ 
tion around an obstacle, refraction in slowly-varying and discontinuous me¬ 
dia, Gaussian beams. 

Minor issues of phase-jumps and corrugations were encountered, and hxed 
respectively by improved selection of zero-search intervals, and by adding 
“ghost cells” outside of the computational domain. 

Future work was outlined and discussed, regarding extensions to anisotropic, 
birefringent, absorbing media such as magnetized plasmas. 
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